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ABSTRACT 


In this work we implement a new method for mapping the Sun’s magnetic field B 
from the photosphere through the corona and interplanetary space. The method 
entails the derivation of B from a scalar potential within a "current-free" annular 
volume bounded inside by the photosphere and outside by a prescribed non- 
spherical "source surface" to which B is made (as nearly as possible) 
perpendicular. As usual we obtain the potential for the part of B that arises from 
currents inside the Sun by fitting an expansion to the observed line-of-sight 
component of B at the photosphere. The new aspect of our work is that we 
introduce a second least-squares fit to obtain the part of B that arises from 
currents outside the source surface. We do this by minimizing the mean-square 
tangential component of B over the source surface. This latter prescription allows 
a nearly arbitrary specification of the source surface, and no particular symmetry 
properties need be invoked for it. We have chosen as our prescription of the 
source surface one that reasonably simulates the expected physical consequences 
of coronal MHD effects. The magnetic field exterior to the source surface is made 
perpendicular to the source surface and is mapped throughout the heliosphere by 
means of a geometrical construction that appropriately generalizes the Parker 
spiral. Detailed comparisons have been made between this new model and 
observed coronal and interplanetary structures. The equatorward inclination of 
coronal helmet streamers is modeled much better with a non-sphericai source 
surface than with a spherical one. There is clear evidence, however, that 
observational data underestimate the strength of photospheric polar magnetic 
fields. Neutral lines on the source surface (/.e., contours on which the normal 
component of B vanishes) generate heliospheric current sheets outside the source 
surface. In the particular case studied, two separate current sheets existed and 
were carried by the solar wind from the source surface to the boundary of the 
heliosphere. The result was a four-sector magnetic configuration near the ecliptic 
at 1AU, and the observations were in good agreement with this prediction of the 


model. 
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1. introduction 


in this paper we deveiop and impiement a significant new tooi for constructing 
models of the coronal and interplanetary magnetic field from observations of the 
photospheric magnetic field. This tool, the non-spherical source surface, represents 
a significant improvement over previous models of solar magnetic field stiucture and 
IS suitable for the study of the interplanetary magnetic field and its solar origin. 

The study of the global magnetic configuration of the solar corona based on 
theoretical extrapolations of photospheric field measurements began with the work 
of Altschuler and Newkirk (1969) and Schatten et al. (1969). Both of these groups 
used the approach of solving Laplace’s equation in an annular volume outside the 
photosphere in terms of a spherical harmonic expansion with coefficients of the 
harmonic functions depending on the measured magnetic field during an entire solar 
rotation. The neglect of coronal currents, necessary to obtain a unique solution in 
simple closed form, has proven acceptable in most cases (Levine and Altschuler, 
1974; Poletto et al., 1975). A harmonic expansion filling all space outside the 
photosphere would result in magnetic field lines that always return to the sun. This 
would not be an acceptable physical approximation to solar conditions because the 
solar wind plasma distorts the outer corona, creating “open" configurations 
consisting of field lines that begin on the sun but do not return. To overcome this 
difficulty, both Altschuler and Newkirk (1969) and Schatten et at. (1969) 
independently introduced a spherical outer boundary concentric with the solar 
surface. The scalar potential between the photosphere and this “source surface" 
was specified mathematically in such a way as to make the coronal field purely 
radial at the source surface itself. The field exterior to the source surface was 
constructed in accordance with the usual Parker (1958) spiral. The mathematical 
implementation of this concept was not difficult. It resulted in a modeling capability 
that accounted at least qualitatively for the effects of the solar wind on coronal 
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magnetic fields. Along with improvements in observational accuracy and resolution, 
improvements in both accuracy and efficiency In the modeling technique based on a 
spherical source surface have been achieved in recent years (Altschuler et a!., 
1976; Adams and Pneuman, 1977; Riesebieter and Neubauer, 1979). 

Applications of potential field models with spherical source surfaces have 
included studies of coronal .structure at eclipse (Schatten, 1968; Altschuler and 
Newkirk, 1969; Smith and Schatten, 1970), coronal currents (Levine and Altschuler, 
1974; Howard et a/., 1980), interplanetary magnetic Helds (Schatten ef a/., 1969; 
Levine et a/., 1977; Levine, 1978; Burlaga et a/., 1978), radio and MHD wave 
propagation in the corona (Uchida et al., 1973; Smerd and Dulk, 1971; Jackson and 
Levine, 1 980), and coronal holes and other sources of open magnetic fields (Levine, 
1977; Svestka et a/., 1977). The success of these applications has been varied. 
The basic structure of coronal holes and the position/polarization of type III radio 
bursts, for example, are very well i’epresented. However, the accuracy of the 
method is clearly less favorable in cases for which the exact structure of the 
outer corona near the source surface is important, as in the study of coronal 
structure at eclipse, the boundaries of coronal holes, and the magnitude of the 
interplanetary magnetic field. These deficiencies are due at least in part to the 
shortcomings of the spherical source surface as a simulation of the MHD effects 
that cause the structure of the outer corona. Many important physical problems, 
such as the possibility of open field lines in active regions, the orientation of 
helmet streamers, the source of the magnetic flux in the ecliptic plane, and the 
topology and evolution of interplanetary current sheets could be investigated with 
more confidence if the structure of the outer corona could be simulated more 
accurately than the spherical source surface allows. 

As a step in this direction, the concept of a non-spherical source surface was 
introduced by Schulz et at. (1978). This yielded a good solution for the case of a 
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magnetic dipole at the center of the sun. Applicability to solar problems was 
demonstrated by showing a close correspondence with the full MHD solution for a 
dipolar magnetic field having boundary conditions of temperature and density 
roughly appropriate to the sun (Pneuman and Kopp, 1971). This correspondence 
opens the possibility for application of the non-spherical source surface to more 
elaborate magnetic field configurations, /.e., to cases for which an MHD solution is 
not feasible. In the present work we extend the concept of a non-spherical source 
surface to the case of an arbitrary photospheric magnetic field and a more flexible 
specification of the source surface. In addition, we show that the field may be 
extended Into Interplanetary space without difficulty, except from regions where 
the source surface is concave. 
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2. Calculation of the Coronal Magnetic Field Inside a 
Non-Spherical Source Surface 

Following the work of Schulz et al. (1978), we will derive solutions for the 
coronal and interplanetary magnetic field by using a non-spherical source surface. 
The goal is to simulate, on a global scale, physical processes whose effects cannot 
be calculated in detail if we wish to retain the great advantage of using measured 
photospheric magnetic fields as the inner boundary condition. As usual, the source 
surface represents the outer boundary of a current-free region. The magnetic field 
inside the source surface is derived from a potential V that is expanded in 
spherical harmonics, as in the earlier work of Schatten et a/. (1969) and Altschuler 
and Newkirk (1969). The magnetic field outside the source surface is based on a 
geometrical construction that properly generalizes the spirahfield concept of Parker 
(1968), and the normal component of B is made continuous at the source surface. 
Because our source surface is non-spherical, the magnetic field will not (in general) 
be radial at the source surface. The field outside the source surface will 
consequently be diverted somewhat toward the main interplanetary neutral sheet 
(Schulz et a/., 1978) and away from the magnetic poles of the sun. Such diversion 
does not occur with a spherical source surface because (in that case) the 6 
component of B is required to vanish at the source surface and beyond. 

As IS usual in source surface calculations, we assume that the corona is 
current-free within the source surface and derive B from a scalar potential of the 
form 



n=0m=0 


P^(®) (g/l* cosm^ + h™ sinm^) 


+ R 


0 


( r \ 

■^1 P™(0)(gJJ*cosm^ + h"sinm^). 


( 1 ? 
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where r, d, and ^ are the usual spherical coordinates and the P^(9) are associated 
Legendre functions. For th|^case of a spherical source surface as outer boundary 
one obtains V=0 at r=Rg by setting 



in (1). Previous authors have sometimes used the symbol r^ to denote the 
constant Rq in the case of a spherical source surface. We use instead the 
equation r=r^(9,^) to describe the configuration of our non-spherical source 
surface, whose mean radius <'r^> is held equal to the constant Rq that appears in 
equation (1). In order to obtain a set of and h^ from observations of the 
photosphere we tentatively adopt a spherical source surface of (constant) radius 
Rq and use equation (2) to eliminate g" and h^ from equation (1). This 
procedure leaves the g^ and h^ as the only remaining parameters to be aoii’^ted. 

Observations provide values of the line of sight field at r=Ro- The 

above harmonic expansion yields an algebraic expression for the line of sight field 
B^=B^ sin(9+/9) -i- B^ cos( 9+/S), where B=-VV is obtained from equation (1) and |3 is 
the average solar B-angle during the time of the magnetic observations (i.e.. |3 is 
the heiiographic latitude of the observer). By requiring that the integral of 
[B/ -(B^sin(9+/3)-*-B^cos(9-«-/}))] over the photosphere be a minimum with 

respect to variation of the coefficients g" and h^, a set of simultaneous linear 
equations for the coefficients can be obtained. If N is the largest principal index 
of the coefficients g™ and h™ to be calculated, there are (N+1)'^ equations and 
(N+1) unknowns, forming a (N+1) x(N+1) matrix problem. However, symmetry 
properties of the harmonic functions decouple the equations into blocks of N-m+1 
simultaneous equations, reducing the largest number of independent equations 
which must be solved at one time to N+1 (Altschuler et a/., 1976). This makes it 
practical to solve for coefficients g^ and h^ with values of N whose corresponding 
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harmonic functions have a structure matching the scaie of synoptic magnetic 
observations ie.g., up to NsQO for synoptic data of approximately 1" resolution). 

For the case of a non-spherical source surface, we recompute the g,!^ and hj!^ 
as an independent set of coefficients whose values depend not only on the 
photospheric field distribution and mean radius of the source surface, but also on 
the shape and orientation of the source surface. The boundary condition at the 
non-spherical source surface should be that the field Is perpendicular to the 
surface. However, for a non-spherical source surface with V given by equation (1) 
this boundary condition cannot be satisfied exactly over the entire surface; it can 
be satisfied only in the least-squares sense. Therefore, a second least-squares fit 
of the g^ and h^ to the outer boundary condition is required. This second least- 
squares fit is the primary difference between the solution for a spherical source 
surface and that for a non-spherical source surface. 

The components of the magnetic field in the region between the solar surface 
and the non-spherical source surface are 


Br = 


2 (g™ cosm^ + h^sinm^) 


n=0m=0 


N n 


2 2n(^) 

n=1m=0 ' 


n-1 


P™(®) (g™ cosm^ + h™ sinm^) 


cm 


(3a) 


N n ,„_.n+2 


= 2* I r J — (g/J'cosm^ + hJJ’sinm^) 

do 


n=0m=0 


N n 


- 2 2 (^) 

n=1m=0^ 


n-1 hdI" 


dP;;'(0) 

d9 


(g^ cosm^ + h™ sinm^) 


(3b) 
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N n 


" " fRol PnW 

= 2 2 I r J m— ^(g;rsinm<^-h;;’cosm<^) 


n=0m=0 


N n 


n — 1 


^ / M" ■ _ 

^ 2 2 Rq m-^|^(g;|’sinm<^-hJ'cosm«^). 

n=1 m=0 ' ’ 

Notice that we have identified a separate value N (#N) for the maximum principal 
Index of the coefficients and h^. The reason for this will be made apparent 
below. 

We let the non-spherical source surface be defined by r^(9,^} with n{d,<f>) the 
unit normal vector field on the source surface and we let RQ=<r^>. We further 
assume that the coefficients and h™ are known. Because we cannot make the 
field be normal to the source surface at each point, we instead minimize the mean- 
square tangential component of B at the source surface. Defining the tangential 
component of the magnetic field at the source surface by 

I ri xB |^ = B^ - (n *B)^, (4) 

we wish to minimize 


o o l*' M 


with respect to each of the g^ and h^. Setting d<r/dg^, =0, we obtain (see 
Appendix for details) the equation 


N n 2n It 
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N n 2ff ff 


2 


{*■''[(" 0 ) O*’”"”*]} (“'^[(r. 


^ 111 •’^sm0 

liJ |n-f| 


Q, P^f(5)cosm'^ 


Ro 


N tf 


^2 2 

Orf'=1nf’=0 


d9d^ 


(6) 


^ ^ ( rf »+1 n* 

S?' J J [ (" •’ [ (^) Prff(»)oo»>'’*] ) {« V [ (r^I P”.'(«)ooam'*] J 


P^W)co«nP*] j pjl'cejcosm’* I 


r^sm6 

In-M 


d9d^ 


^^2 2 

0 rf= 1 mr =0 


N n' 


^ / rf+1 n' 

wf J I [(«•'' [(-^) Prf!'(«»ln"f*]}(« ''[('^) PjWcosm'*]] 





n’ 

P^t(0)cosm'^ 


p 

r'^sinO 

d0d^ 

I yv ^ I 

nr] 


A similar expression is obtained from the condition d<r/dh^ =0, except that the 
cosm'^ terms in equation (6) are replaced by sinm'^. 

Together, equation (6) and its analog represent N(N+2) simultaneous linear 
equations for the unknowns gjIJ’ and h™. There is one equation for each n' for 
m’ = 0, and there are two equations for each combination of indices 
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(n*; m')=(1, 2 N; 1, 2 n' ). The coefRcients of the unknowns g™ and h^ 

are the integrals on the left hand side of equation (6) and its analog. These form 
a matrix with N (N+2) x N (N-t-2) elements (see Appendix for details). The right hand 
side of each equation consists of the integrals on the right hand side of equation 
(6) and its analog. These form a column vector of N(N+2) elements. Because any 
continuous with n^=ri*r^O everywhere can be specified as the source 

surface, it is clear that no general symmetry exists in equation (6) which reduces 
the number of independent equations in all cases. The harmonic functions are not 
mutually orthogonal over a non-spherical source surface because both r and n -f 
depend on d and ^ in the general case. There is not even enough symmetry in 
general to decouple the linear equations corresponding to distinct values of n or m. 
Thus the full N(N-i-2)xN(N+2) matrix must be calculated and solved to obtain a 
solution for a set of g^ and h^, and it is computationally prohibitive to set N=N 
We can possibly afford to take N^go in deducing the g^ and h^ from line-of-sight 
data, but we cannot afford to take N 1 5 in deducing the g^ and appropriate 
to a non-spherical source surface. The former operation requires us to compute ~ 
N /3 matrix elements (/.e., 2.4x10 elements for N=90), while the latter 

operation requires us to compute matrix elements (/.e., 4x10^ elements for 
N=14 and 6.4x10^ elements for N-90). 

As a compromise between thrift and resolution, we typically take N=25 (with an 
integration grid of 60 latitude zones by 120 iongitude sectors) and N=6 with a grid 
of 18 X 36). Further discussions of this solution are given in the Appendix. 
Calculating the matrix elements (which requires the most computation time) and 
solving the system of linear equations (which uses less than 6% of the computation 
time) typicaliy takes 20 times longer for the g™ and h^ up to N=6 than for the g^ 
and h^ up to N=25 (even with the finer integration grid used at the photosphere). 
This ratio is insensitive to the computing power applied, being essentially invariant 
between a POP 1 1/60 minicomputer and a CRAY-1. Calculation of such a set of g^ 
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and can be completed in approximately one hour on a PDP 1 1 /60 or 7 seconds 

-4 

on a CRAY-1. This time Increases approximately as 2N . For example, a solution 
for N=14 (224 x 224 matrix) requires 7 minutes on a CRAY-1. 

Because the contribution of the and h^ terms in equations (3) is non-zero 
at r=RQ, the full solution derived above will not satisfy quite the same 
photospheric boundary condition that the solution for a spherical source surface 
satisfies. A fully self-consistent model that satisfies both boundary conditions at 
once would seem to require a simultaneous least-squares solution for al! the 
coefficients when the source surface is non-spherical. In the present application of 
our model, however, the configuration of the source surface Itself depends on the 
and h^ (see below). Thus, any simultaneous fit for the g^, h^, g^, and h^ 
would be highly nonlinear in the g^ and h^. The implementation of such a scheme 
would be computationally prohibitive. However, we do not expect the deviations of 
q™ find h^ from (2) to be so targe as to invalidate our determination 
(described above) of the g^ and h^ from line-of-sight data at the photosphere. 
Our confidence rests on the fact that the second line of equation (1) contributes 
little to the field at r=RQ in any event, /.e., about one part in (Rq/Rq^*^ ' 

harmonic number n. Thus, in the numerical implementation of equation (6) and its 
analog, the coefficients g^ and h^ which are used are those defined by the 
measured photospheric field distribution but with a spherical source surface at 
r=<rg>. This use of a spherical source surface to determine the "photospheric" 
coefficients g^ and h^ is easily implemented, since it modifies but slightly the 
algebraic form of B/ from that used by Altschuler et al. (1976). The values of g^ 
and h|!^ thus determined are also used in the calculation of magnetic field 
components by means of (3). 
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3. Specification of the Non-Spherical Source Surface 


The formulation of Section 2 is applicable to any non-spherical source surface, 
subject only to the restrictions that r^(9,^) is continuous and that n^=n-r:^0 
e«rerywhere. For practical application it is necessary to choose a definition of the 
non-spherical source surface that is both physically motivated and computationally 
acceptable, in order to distinguish clearly between the physical prescription of the 
source surface and its application to the solution for B, we define a general 
analytic function F(r,9,^) and require the source surface to be a surface of 
constant F. We have chosen an initial definition of F which is based on physical 
expectations and which has particular computational advantages. We expect that 
refinements will be made in the physical prescription used to choose a source 
surface as the usefulness of this method is explored. 

The purpose of a non-spherical source surface is to simulate the Implications of 
having a non-radial solar wind at and beyond the outer boundary of the free 

region. The non-radial solar wind is a fully expected consequence of MHD 
(Pneuman and Kopp, 1971). Schulz et a/. (1978) suggested that the magnetic field 



n=1 m=0 


(g™ cosm^ + h™ sinm^) P^(^) 


(7) 


derived solely from internal sources, /.e., from the first line of (1), be considered as 
the basis for specifying the source surface. By choosing F=B , Schulz et a/. 
(1978) specified the source surface as an isogauss of (7), /.e., as a surface of 
constant | B | , and thereby obtained a rather good simulation of the Pneuman-Kopp 
MHD solution for the case of an internal solar dipole. 

For more realistic cases, however, we have found too often that an isogauss 
of (7) is not a convex surface. In other words, it happens too often that the 
source surface specified by Schulz et a/. (1978) is indented in places. This 



condition seriously complicates the geometrical construction of B outside the 
source surface, because it permits solar-wind streams from different parts of the 
source surface to intersect. Intersecting streams may well be present in the real 
solar wind, but they impose mathematical difficulties that we wish to avoid m the 
present model. Accordingly, we have re-opened the question of how to specify an 
optimal source surface and have discovered a suitable compromise between the 
spherical source surface and the isogauss of | B | . 

Our new procedure is based on a more general specification of F, namely 


F =aB^ + (1-o)Bf 


<r?> 


( 8 ) 


in which a and v are parameters to be chosen. We have fixed for 

convenience The effect of choosing the parameter a<1 is to generate a surface 
of constant F v'lat is smoother than an isogauss of B, but has the same average 
radius. Thus the case a=0 produces a spherical source surface of radius Rq, while 
the case a=1 produces a source surface that is an isogauss of B with average 
radius <r^>=RQ. The constants B^ and <r^> are determined in an initial step by 
setting a=1 in (8) and finding the values of B and <’r^> that correspond to an 
isogauss of mean radius <r^>=RQ. Thereafter we generate surfaces of constant F 
for any desired values of a<1 by holding <r^>=RQ, but we do not recalculate the 

~2 V 

constants B^ and <r^> 

The outward normal to the source surface is the unit vector in the direction of 
-VF, /.e., the unit vector n =-VF / 1 VF | . Because F is expressed analytically by 
equation (8) with the help of equation (7), we can also express n analytically, and 
there is no need for numerical differentiation to calculate ri. 

Some implications of selecting a value of a other than asQ (Altschuler and 
Newkirk, 1969) or o=1 (Schulz ef a/., 1978) are illustrated in Figures 1-4. Figure 1 
shows a longitudinal section of the source surface corresponding to a photospheric 
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Figure 1 Longitudinal section of the source surface corresponding to a 
photospheric dipole for three different values of a. The average 
radius of each surface is 2.68 Rq. 
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Figure 2 Longitudinal sections of the source surface calculated from the 
measured photospheric magnetic field for Carrington rotation 1602 for 
three different values of a. Where the a=0.6 surface is not shown, it 
essentially coincides with the a=1 surface. The average radius of 
each surface is 2.3 Rq. 
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Figure 3 Comparison of magnetic field lines calculated by the present model 
(solid lines) with field lines calculated by the MHD model (dashed 
lines) of Pneuman and Kopp (1971) for different values of a. The 
footpoints of our computed field lines were chosen to match those of 
the Pneuman-Kopp field lines. The location of the source surface is 
indicated by dotted lines. 
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Figure 4 Contours of constant r^ for Carrington rotation 1 602, with two 

w 

different values of a, in synoptic chart format. Contour lines are 
shown at intervals of 0.1 Rq; the contour at rg=<rg>=2.3 Rq is 
thickened. The arrows are projections of the unit normal vector onto 
the unit sphere. 
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dipole (g^^O) for each of several values of a. The value of <r^> is held 
approximately equal to 2.58 Rq In each case. Figure 2 shows certain longitudinal 
sections of the source surface corresponding to the g^ and hjl^ deduced from the 
measured line-of-sight field at the photosphere for Carrington rotation 160?. The 
values of a are the same in Figure 2 as in Figure 1, but the value of \r^> is held 
equal to 2.3 Rq Figure 2. We see that a choice of a somewhat less than unity 
can often eliminate concave areas (indentations) that occur on the source surface 
for a=1. 

For the present, we assume that the optimal value of a in general can be 
determined by reexamining the case of a photospheric dipole. Thus, we have 
constructed a source surface of mean radius <’r^>=2.58RQ for each of several 
values of a and have carried out the other steps described above for the case in 
which g^ is the only nonvanishing coefficient among the g^ and h^ All the 
coefficients g^ and h]![' that should vanish in this case (/.e., all except and 

'ig for N=6) indeed turned out to be very small. We computed B from (3) and 
plotted certain of the field lines on a meridional plane. The results are shown in 
Figure 3 as solid curves. The dashed curves in Figure 3 are the field lines 
obtained by Pneuman and Kopp (1971) in their MHD study of the photusphenc- 
dipole case. Our field lines were chosen to have the same heliographic latitudes at 
r=RQ as the Pneuman-Kopp field lines. 

The choice a=1 in Rgure 3 corresponds to the solution of Schulz Ht af (1978), 
although our Figure 3d differs somewhat from their Figure 2. Such minor differences 
should be expected, since (a) our <r^>=2.58 Rq exceeds their <r^ >=2.547 Rq by 
about 1%, (b) our integration grid in equation (6) is necessarily coarser than theirs, 
(c) our N=6 IS necessarily smaller their N=1 5, and (d) our procedure for tracing field 
lines differs somewhat from theirs. Regardless of these differences, it is ciear that 
the choice of «=1 implicit in the work of Schulz et a/. (1978) produced a good 
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simulation of the MHD results of Pneuman and Kopp (1971), whereas the 
adjustability of a in the present work offers the prospect of Improved simulation. It 
IS particularly important in this context to model well the asymptotic directions of 
the field lines, because one of our mam motives for developing the nomspherica' 
source surface has been to extend the model as accurately as possible throughout 
the heliosphere. It is clear from Figure 3 that the choice «=1 provides too much of 
an equatorial deflection of field lines from the radial direction outside the source 
surface, whereas the choice a=0.4 provides too little. The four panels of Figure 3 
all show rather good agreement between our model and the MHD results Inside the 
source surface (dotted curve) in each case. In selecting an optimal a, we do not 
rely heavily on comparative details in the vicinity of the cusped MHD field line. Our 
model Is current-free within the source surface and so inherently produces a Y*type 
neutral line, whereas MHD produces a cusp-type neutral line (Sturrock and Smith, 
1968; Pneuman and Kopp, 1971). Taking these various factors into account, we 
conclude from Figure 3 that a choice of a in the range 0.6-0.8 best simulates the 
MHD solution overall. We have seen from Figure 2 that this range of a tends to 
yield a convex source surface when used with photospheric data, as desired. 

Although we will restrict our further examples to isogauss surfaces with a In 
the range 0.6-0.8 (cf. equation (8)), it should be remembered that this choice is 
based on comparison with an MHD solution that is azimuthally symmetric and 
represents an isothermal corona. Both of these restrictions are serious 
compromises with the known physics of the solar corona. Thus, It would not be 
surprising if further progress in simplified MHD models of the corona pointed the 
way to different ranges of a, or, more likely, to different ways of specifying a 
three-dimensional non-spherical source surface. 

A global view of two of the source surfaces profiled in Figure 2 is given in 
synoptic chart format in Figure 4. Here the contours of constant heliocentric 
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distance on the source surface are incremented in steps of 0.1 Rq, and the 

contour r=<r^ >=2.3 Rq is heavier than the others. The arrows represent 

projections of the unit normal vector n on the unit sphere, and so have lengths 
p pH 

proportional to (n^ +n^) . The shorter projections thus correspond to the more 
nearly radial orientations of the unit normal. Concave areas are not readily 
Identified in this format, but the presence of fewer contours for oe=0.6 suggests 
that the source surface is more nearly spherical (as expected) for oe=0.6 than for 
a=1. Centers of convergence (toward which the arrows point) do not necessarily 
indicate concave areas but do indicate regions above which |B| diminishes less 
rapidly than Parker (1958) suggests with heliocentric distance beyond the source 
surface. Centers from which arrows emanate identify regions above which |B[ 
diminishes more rapidly with r than Parker (1958) suggests. The variation of B 
interior to the source surface enables open field lines to map the entire source 
surface onto a rather small fraction of the photosphere (Levine, 1977). Thus, the 
variation of |B| with r inside the corona must be considerabiy stronger than 
inverse-square on the average. Most coronai magnetic modeis possess this iast 
property. Ours offers the further prospect of modeiing the variation of ]B|, and thus 
of flux-tube area as a function of r (Nemey and Suess, 1975; Suess et al., 1977; 
Levine, 1978), ail the way from the photosphere to the outer boundary of the 
heliosphere. 


21 



This Page Intentionally Left Blank 



4. Coronal Magnetic Field Models inside the Non-Spherical Source Surface 
4.1 Accuracy of the Solution for and 

In Section 2 it was pointed out that the least-squares solution for the 
coefficients and h^ for a non-spherical source surface will not, in general, 
satisfy the boundary condition BxnsO exactly at all points on the source surface. 
In addition to accuracy in the above sense, we consider in this section the 
precision of the numerical procedures used to solve equation (6) for the 
coefficients g^ and h^. (The stability of the solution Is discussed in the 

Appendix.) 

The precision of the numerical solution of equation (6) can be checked by 
calculating coefficients in cases for which the answer is known by some 

independent means. One such case is a spherical source surface, which can be 
considered a degenerate case of equation (8) with a=0. In this case the 

coefficients g^ and are given exactly by equation (2). Our procedure entails 
the numerical computation of matrix elements over a finite integration grid, and so 
our results for g”^ and h^ do not agree exactly with (2). However, calculation of 
g^ and h^ using equation (6) with a source surface given by 

a=0(/.e., r^(9,i^)=RQ; n=r) satisfies equation (2) to within 0.1% in aii cases that 
we have examined. 

A more stringent test is the case of a photospheric dipole with a source 
surface that is an isogauss of B (/.e., a=1). Schulz ef aJ. (1978) calculated the 
coefficients g® (n=1, 3, 6, . . ., N) for this case by obtaining a solution that 
assumed azimuthal and north-south symmetry and that specified r^(9,^) and n 
analytically. They obtained accurate numerical results for g^/g^ (for various 
values of N), and needed to solve at most (N+1)/2 simultaneous linear equations 
for this case, because the symmetry of the problem required that g™ and h^ 
vanish unless m=0 and n is odd. Our present code does not constrain any of the 
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coefRcients or in (1) to vanish. Numerical solution of equation ( 6 ) under the 
conditions discussed in Section 2 (r^ defined on a grid of 18x36 values; N= 6 ) 
gives agreement with the values of g^, g^, and gg calculated by Schulz et a/. 
(1978) to within 0.30% for g®, 2.60% for gg, and 3.74% for gg. The other 
coefficients should be zero for this dipole case. For m^O this is true in our 
numerical calculation by at least six orders of magnitude compared to g^. For m=0 
the average of the absolute values of the coefficients g^i 94 ) and gg is 0.56% of 
the average of the absolute values of the the three non-zero coefficients. 

We have investigated separately the effects of increasing N and cf doubling 
the number of source surface grid points in our numerical calculation. Each of 
these can improve the agreement with the results of Schulz et eJ. (1978), but 
there appears to be a limit to the agreement of the n=1, 3, and 5 coefficients at 
about one-haif to one-third of the percentages given above. Moreover, with Ni.14, 
the largest value tested in our numerical scheme, the agreement was not 
satisfactory for odd values of n^9. We conclude that, with the present numerical 
procedure, optimal solutions (In terms of overall agreement with the results of 
Schulz et al., 1978) are obtained for N =6 or 8 and an integration grid at least as 
fine as 18 X 36, and preferably 36 x 72. The increasing size of the fairly sparse 
matrix for larger N (e.g., 224 x 224 for N=14) seems to introduce numerical 
difficulties that adversely affect the convergence of the coefficients having higher 
values of n. 

The minimization of <f for a non-spherical source surface does not (in general) 
cause <r (equation ( 6 )) to vanish exactly, but it should cause the strong inequality 
|B X n « (Rq/Rq)® [(g ® )^ + (g .j )^ + (h .J )^] to hold over most of the source 
surface. Thus, we expect to find v « Air (Rq^/Rq^^ 1 1 1 

the optimal expansion coefficients g^ and h^, and we expect to find the integrand 
I B X n I to be distributed like random noise over the source surface. Our 
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numerical results largely confirm these expectations. For the case of a 
photospheric dipole with or=0.8 (see Figure 3c) we found 

<r=0.0175x4ir(Ro®/Ro^^ I B x n || < 0.087 (Rq/Rq)® <9 ? 

everywhere. However, because this special case was azimuthally symmetric, we 

o 

found no variation of |B xn|^ with longitude. Instead we found, as did Schulz ei 

p 

a/. (1978), that contours of constant |Bxnjg divided the source surface into 
latitudinal bands. 

In more realistic cases, e.g., for the and h^ deduced from Ime-of-sight 
photospheric field data for Carrington rotation 1602, the azimuthal asyrmuetry of 
the source surface (see Figure 4, above) causes the distribution of | B x n | ^ to 
appear essentiaiiy disorganized with respect to 6 and 0, as desired. Moreover, the 
vaiues of <r turn out to be small, satisfying the inequalities proposed above. We 
have verified that the above inequalities are satisfied more strongly in all cases 
when the number of grid points on the source surface is doubled and (with 
°r>mewhat less improvement) when larger N is used. Thus, we have good 
confidence in the utility of our model for application to actual coronal and 
photospheric data, as well as to idealizations of the sun’s magnetic field. 

4.2 Comparison with Coronal Data 

We have chosen the Skylab ATM era as the best time period for comparison 
with coronal data because of the completeness of the data set. Because the most 
visible effect of incorporating a non-spherical (rather than spherical) source surface 
into coronal magnetic models is that the calculated helmet streamers now can be 
non-radiai, we have selected the ATM coronagraph data as the most relevant ATM 
data for comparison with the model. In particular, the total eclipse of 30 June 
1973 yields the best data for comparison because there are also ground-based 
photographs that show coronal structure down to r=RQ. Furthermore, there was 
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relatively little evolution of the sun’s magnetic field during this time period (Levine, 
1 978). 

The non-spherical source surface adopted for the eclipse time period is shown 
in Figure 5 in the same format as Figure 4. For the calculation of this surface and 
the coefRcients g™, h™, g™ and h^ photospheric magnetic data from Carrington 
longitude 180° of rotation 1602 to longitude 180° of rotation 1603 were used. In 
calculating the non-spherical source surface we have been guided by the 
experience of the dipolar case (see Figure 3 above) to choose a value of a=0.8 in 
equation (8). The value of <r 3 > was chosen to produce a source surface height 
that matched the tops of the helmet streamers visible at the 1973 eclipse. The 
average radius of the source surface In Figure 6 is 2.72 Rq and the extrema are 
1.92 Rq and 3.17 Rq. 

The effect of making the source surface non-spherical can be seen by 
comoaring Figure 6, which is the ground-based eclipse photograph, with Figure 7. 
Figure 7a is based on a potential field model with a spherical source surface at 
r=2.72 Rq, which is the average radius of the non-spherical source surface used 
to construct Figures 6 and 7b. In Figures 7 and 8 starting points for tracing field 
lines have been chosen to emphasize the boundary field lines between open and 
closed structures. These field lines should lie immediately over helmet streamers 
and on the outer boundaries of coronal holes or other coronal structures. This 
choice is accomplished by starting to trace from a small distance on either side of 
polarity reversal contours of vanishing B • n on the source surface. Moreover, in 
order to emphasize the direction of the magnetic field lines beyond the source 
surface, each field line has been traced for 0.6 Rq beyond the source surface. 
For a spherical source surface this extension is radial, but for non-spherical source 
surfaces it is along the local normal vector n. Superimposed on these field lines 
are the field lines that start at the 400 locations on a grid of 1 8 x 36 points at 


26 



LATITUDE 



0 ° 90 ° 180 ° 270 ° 360 ° 

CARRINGTON LONGITUDE 

Figure 5 The source surface and unit normai vector for Carrington rotation 
1602/3 (16 June 1973 to 12 Juiy 1973). The format is the same as 
in Figure 4. The thickened contour ilne is rg=2.72RQ (<’rg>=2.7RQ 
for this case). 
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Top: a=0, /.e., the model with a spherical source surface. Bottom: 


«s=0.8, /.e., the model with a non-spherical source surface. Dashed 

fiAlri ItnAici lif> hf^hinri thei nlane nf the . <kv. See tevt for de.«5r;ribtinn of 


the algorithm for choosing which field lines to plot. North to the top, 
west to the right. 
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r=Ro which have the largest values of Bj.. These field lines emphasize the closed 
structures that underlie the helmet streamers. 

As expected, the primary difference between the two models is in the direction 
of the helmet streamers. Only a subtle difference in the corona, this will result in 
major differences in the predicted magnetic field strength throughout the 
heliosphere when the models are extended to r^lAU. A total of fifteen features 
have been identified extending beyond 2.5 Rq on the eclipse photograph of Figure 
6, and on the corresponding ATM coronagraph photograph (R. H. Munro, private 
communication, 1980). Features that extend radially outward from the spherical 
source surface (Figure 7a) correspond to similar features that extend normally 
outward from the non*sphericai source surface (Figure 7b). For each feature the 
deviation from the radial direction was measured on the photograph and calculated 
from the direction of n at the source surface. The average equatorward deflection 
of the observed coronal streamers was 6.2°, with a rather large dispersion of 5.9°. 
The model streamers for a non-spherical source surface showed an equatorward 
deflection of 2.1° with a dispersion of 3.2°, indicating approximately a 30% 
improvement in this important aspect of coronal field modeling. Histograms of the 
equatorward deflections of the chosen features are shown in Figure 9. 

In addition to the remaining 4.1° of equatorward deflection that the non- 
spherical source surface did not model, there are still other major discrepancies 
between the model and the observed coronal structure. The most glaring 
discrepancy is that the model has placed the large streamer in the southwest 
quadrant in the wrong location. It is calculated to be about 30° poleward of where 
It should be. This mislocation is far too great to be explained by a poor choice of 
the model parameters or by other inaccuracies of the model. Furthermore, the 
mislocation extends down to the very surface of the Sun, indicating that the 
source of the discrepancy lies with the input data. Because this particular 
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EQUATORWARD DEFLECTION (OBSERVED MINUS CALCULATED! 


Figure O Histograms of the differences (observed minus caiculated) In the 
equatorward deflection of 16 coronal features observed outside the 
source surface: (a) Spherical source surface, /.e., a=0; (b) Non- 
spherical source surface, asO.8; (c) Non-spherical source surface, 
a=0.8, with added polar flux. Because the calculated equatorward 
deflection for the spherical source surface is zero, panel (a) Is In fact 
a histogram of the observed deflections. 
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streamer overlies a polar crown prominence, we suggest that the polar field 
strength has been underestimated by the magnetograph observations. 

This is not a new hypothesis. It has been suggested by Pneuman et al. 
(1978) on the basis of similar data and reasoning. It has also been suggested by 
Suess et al. (1977) on the basis of coronal hole modeling, and by Svalgaard et at. 
(1978) on the basis of direct magnetograph data. The difficulties that 
magnetographs have when measuring polar fields has been a topic of extensive 
discussion for many years (cf. Howard, 1977, and references therein). In this 
latest case, we seem to have clear evidence (/.e., the mislocation of a polar crown 
prominence) that the polar fields are indeed underestimated. 

To test the hypothesis for this particular case, we have added polar magnetic 
flux to the magnetograph data and recomputed a new model (/.e., a new source 
surface with the same es before and new coefficients g^, h^, g^, h^). 

Polar flux has been added in accordance with the model proposed by Svalgaard et 
a/. (1978). Specifically, poleward of 70° latitude we have added to the 

O 

magnetograph data the radial magnetic field AB=±BQr cos 9, with Bq= 1 1.56 G, the 
sign being chosen so as to reinforce the predominant polarities of the respective 
polar caps. We do not necessarily endorse this particular model nor imply that it is 
the only model that yields an acceptable fit with the coronal data. We have 
selected it because it is complete, recent, and based on magnetograph data only. 
There is probably a range of models for this "missing polar flux" that would produce 
similarly acceptable results. 

Field line plots for the model with the added polar flux are shown in Figure 8. 
Figure 8a is based on a spherical source surface and Figure 8b on a non-spherical 
source surface (a=0.8) with the same added polar flux. As expected, the added 
polar flux has relocated and reoriented the helmet streamer m the southwest 
quadrant, in much better agreement with the observed streamer. Also, the average 


33 



equatorward deflection of all the chosen features Is modeled much better with the 
added polar flux. The difference between observed and computed deflection is 
reduced from 4.1° without the added polar flux to 1.5° with added polar flux. The 
agreement is still not as good as one would like, but the present comparisons are 
sufficient to show that the addition of polar flux is necessary and sufficient to 
model polar crown streamers reasonably. 

The large north polar coronal hole visible in Figures 6, 7, and 8 is the same 
hole that was studied extensively by Munro and Jackson (1977). A comparison 
can be made between the coronal hole boundaries measured by Munro and Jackson 
from coronagraph data, x-ray photographs, and K-coronameter data and the 
boundaries that would be inferred from the present models. This comparison Is 
shown in Figure 10, where we have used the model that contains the added polar 
flux, /.a., the model of Figure 8b. This model fits the observed boundary extremely 
well on the east limb, but the fit is poor on the west limb. The position of the 
predicted boundary is altered only slightly by making the source surface non* 
spherical rather than spherical, but the added polar flux improves the comparison 
greatly. 

There are at least three factors which could contribute to an explanation of 
the discrepancy at the west limb. The first is the possibility that the photospheric 
magnetic field evolued signiflcantiy in the region underlying the west limb during the 
period between its central meridian passage and the eclipse (approximately one 
week). The second is the possibility that there may be closed magnetic structures 
in the corona which do not contain enough plasma to emit detectable amounts of 
visible or x-radiation. An example of such a situation would be the northern portion 
of the helmet structure on the northwest limb in Figure 8b. The third Is the 
possibility that the polar regions may contain even more unobserved magnetic flux 
than proposed by Svalgaard e< a/. (1978). (The total additional polar cap flux in 
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Figure 10 Comparison of the boundary of the north polar coronal hole as 
predicted by the present model (solid lines) and as measured (dashed 
lines) by Munro and Jackson (1977). 
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their model amounts to 3.24x10 Mx.) The continued dilTerence between the 
modeled and the observed streamer In the southwest quadrant argues in favor of 
this explanation. 

An additional factor in the comparison of the model with coronal observations 
indicates an inherent limit to the applicability of potential field models themselves. 
The models of Figure 8 consistently calculate helmet streamers that are broader 
than observed. This is true In both equatorial and polar regions. A subsidiary 
result is that the calculated areas corresponding to the footpoints of open field 
lines are consistently smaller than the corresponding areas of observed coronal 
holes. This difference was studied by Levine (1977) using a spherical source 
surface, and it persists in the present work with non-spherical source surfaces. 
Because this difference is found at all latitudes, it cannot be ascribed to 
unobserved polar flux. It seems more likely that the observed helmet streameis 
are narrower than the calculated ones because there are current sheets at the 
boundaries between open and closed field lines, even Inside the source surface. 
Such current sheets are, of course, not included in potential field models. One 
effect of such current sheets could be to make helmet streamers narrower than if 
the currents were not present. Physical reasons for expecting such current 
sheets are discussed in Section 6. 
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5. Magnetic Field Outside the Non-Spherical Source Surface 


Our inodel magnetic field outside the source surface is given by a geometrical 
construction. In this construction we differ from Parker (1958) in detail, but not m 
spirit. Indeed, our basic principle is that solar wind plasma elements follow 
rectilinear trajectories upon emerging from the source surface. However, the 
trajectories are not radial because (a) the source surface is not spherical and (b) 
the source surface rotates with the sun. The geometrical construction requires us 
to transform between the non-rotating frame (in which the trajectories are 
rectilinear outside the source surface) and the co-rotating frame (m wh'c^ the 
plasma elements flow parallel to B inside the source surface. 

In the co-rotating frame the solar^wind velocity v Is ideally normal to the source 
surface, as is B. (Within the source surface we have implemented this ideal only in 
the least-squares sense, but the procedure does not invalidate the ideal.) It 
follows th 2 t v=” Vjj. We assume for simplicity that v^ is constant over the source 
surface and adopt the value v^=320 km/sec recommended by Hundhausen (1970). 
(One can permit v^ to be a function of latitude, but not longitude, if the source 
surface is azimuthally symmetric. In general, the permitted variation would be more 
subtle than this, and beyond the scope of the present work to describe. 
Interplanetary shocks result from the failure to impose appropriately restrictive 
conditions on the variation of v^ over the source surface.) It follows that the solar 
wind velocity in the non-rotating frame is given by 

u = n v^ + (jgx r)g , (9) 

where Q is the angular velocity of the sun (corresponding to a siderial period 
2 it/Q« 25.1 41 4 days) and the subscript s denotes evaluation at the source 
surface. 

It proves convenient to work in rectangular coordinates and to describe the 
trajectories parametrically by means of the time t after emission of the plasma 
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element from the source surface. We let the co-rotating system (x, y, z) coincide 
with the non-rotating system (x*, y', z') at time t^O (see Figure 11). The velocity 
of the plasma element In the non-rotating system Is given by (0) as 

u = (njj v„ - Qy^) x' + (ny v„ + Ox^) y’ + n^ v,, z' (1 0) 

at time t^O. We assume that the velocity remains constant for a|| t along the 
rectilinear trajectory. We assume further that 6 is a constant vector parallel to z' 
/.e., that differential rotation does not apply to large-scale magnetic fields. It 
follows from equation (10) that the trajectory of a plasma element is given 
parametrically by 


x'-X 3 + (nx v„ -QygJt 
y‘ = y8+(ny v^+OxgJt 

z'-Zj+nzVnt 

for t>0, where n^=(il * ^l),. However, it follows from Figure 1 1 that 
x* = X cos Ot - ysin Ot 


y* = X sin Ot + ycos Ot 


(lla) 

(llb) 

(llc) 

(12a) 

(12b) 


and z'sz. Substitution of (12) in (11) yields a pair of simultaneous linear equations 
for X and y, and it follows that 

X = [Xg + (n^ Vj, - Qyg)t] cosOt + [y^ + (Py v,, + Qxg)t]slnOt (13a) 


y = [Vs + (Hy v„ + Oxg)t] cosOt - [X3 + (n,j v„ - Qys)t] slnOt (13b) 


along a plasma streamline in the co-rotating frame. Moreover, it follows from 
Parker's principle that (13) also constitutes the parametric equation of a field line 
outside the source surface. 
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Figure 11 



Relation between the rotating coordinate system (x. y, z) and the 
non-rotating coordinate system (x’, y', z’). The z axis coincides with 
the z' axis and extends out of the plane of the paper. 
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While equation (13) can be used to trace aq^ tfeld line from the source 
surface to the outer boundary of the heliosphere, we use (13) in the present work 
to map the interplanetary current sheets associated with neutral lines on the 
source surface during June-July 1973. We use the photospheric data from 
rotations 1602 (last half) and 1603 (first half) with the added polar flux (as in 
Figure 8). There were two neutral lines (/.e., contours of vanishing B • n ) on the 
source surface during this (Skylab-ecllpse) epoch. These are shown In synoptic- 
chart format in Figure 12a as solid curves. By means of equation (13) we have 
mapped these neutral lines from the source surface to heliocentric spheres of 
increasing radius. The loci of the current sheets at r^lORQ and r^SO Rq are 
shown in Hgure 12a as dashed and dotted curves, respectively. Figure 12b shows 
the locus of each current sheet at r-IAU (solid curves) and r=2AU (dashed 
curves). The displacement of these loci in longitude at interplanetary distances is 
" c ''rs&auc;nc& of the spiral configuration of the magnetic field. The latitudinal 
displacement of current-sheet loci in the Inner heliosphere (r^lORo) is a 
consequence of significantly non-radial flow of the solar wind there. 

The polarity of the Interplanetary magnetic field observed at r=1AU (King, 
1976) is indicated on Figure 12b. It can be seen that two of the "sector boundary 
crossings" are very well modeled by our projection of the current sheets to 1AU, 
whereas the other two crossings are not well modeled. These latter crossings, 
however, are known to have enclosed a very high speed solar wind stream. We 
can see from this the need for properly modeling the variation of v^ over the 
source surface (rather than taking v^ as constant, as we have done for 
convenience in the present work). Moreover, the principal effect of a non-spherical 
source surface is on the latitudinal structure of the current sheet and comparison 
with observed ecliptic polarities is a limited test of the capabilities of the method. 


' 9 ' 
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Figure 12 Location of the interpianetary current sheets for Carrington rotation 
1602/3 at different heliospheric distances in synoptic chart format. 
Panel (a): r<100RQ. The latitudinal displacement of the current 
sheet results from the non-sphericity of the source surface and takes 
place mostly within 10 Rq. Panel (b): r>100RQ. The longitudinal 
displacement perceptible in both panels results from the spiral 
configuration of the magnetic field. The measured polarity of the 
magnetic field at r=1AU in the ecliptic (King, 1975) is indicated in 
panel (b). This polarity should change sign upon traversal of the solid 
line. 
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6. Summary and Discussion 


We have developed the concept of the non-spherical source surface into a 
well-defined operational method for simulating the magnetic structure of the real 
solar corona and the real heliosphere. The method follows precedent by deriving 
the coronal B field inside the source surface from a scalar potential and the 
heliospheric B field outside the source surface from a geometrical construction. 
However, the non-spherIcal character of the source surface enables our method to 
simulate the non-radial orientation of helmet streamers and other characteristics of 
solar wind flow. We know from past analyses (Schulz, 1973; Schulz et at., 1978) 
that some such recognition of the non-radial component of solar wind flow in the 
inner heliosphere is necessary in order to account for the observed magnitude of 
|B| in interplanetary space. 

Our source surface is specified by equation (8) and contains an adjustable 
parameter a such that 0 ^ a ^ 1 . The limit a=0 corresponds to the spherical source 
surface of Altschuler and Newkirk (1969) and Schatten et at. (1969); the limit <t=1 
corresponds to the non-spherical source surface of Schulz et at. (1978). The range 
0.6 ^ a ^ 0.8 yields optimal agreement between our model and the full MHO solution 
of Pneuman and Kopp (1971) for the case of a photospheric dipole. It is our 
working hypothesis that the same range of a (i.e., a»:0.6-0.8) is appropriate to the 
study of the real corona also. 

Using photospheric magnetograph data centered on the Skylab-era eclipse of 
30 June 1973, we have found that our present method adequately models the non- 
radial orientation of helmet streamers in eclipse and coronagraph photographs if 
additional magnetic flux is superimposed on the polar regions. However, observed 
helmet streamers are consistently narrower in latitude than those calculated by our 
method, and observed coronal holes are often larger in area (cf. Levine, 1977) 
than the footpoint regions of open field lines calculated by our method. These 
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discrepancies are tentatively attributed to the neglect of coronal current sheets 
that are likely to occur (even at r<r^) at the boundary between closed and open 
Held lines in consequence of the discontinuity in plasma pressure, it Is known, 
after all, that coronal holes correspond to regions of lower plasma density and 
temperature than the coronal arcades that form helmet streamers (e.g., Bohlin, 
1977). Thus, the discrepancies between observed and calculated coronal hole 
areas and helmet streamer widths may well be intrinsic to our model, as well as to 
other models in which B is derived from a scalar potential within the source 
surface. 

Although mathematical dlfltculties prevent us from enforcing the desired 
boundary condition |B xn IbO at every point on the source surface, our procedure 

O 

of minimizing the mean value of | B x n | proves to be quite successful. Field-line 
plots (e.y.. Figure 3) show that 6 is nearly perpendicular to the source surface 
everywhere except near the neutral lines from which Interplanetary current sheets 
(e.g., Figure 12) originate. However, the derivation of B from a scalar potential in 
the corona necessarily causes neutral lines on the source surface to be Y-type 
rather than cusp-type (Sturrock and Smith, 1968; Pneuman and Kopp, 1971). Thus, 
one could not expect the last closed coronal field line to be perpendicular to the 
source surface there, even to a spherical source surface. Our field-line plots 
(Figures 3, 7, and 8) show that the least-squares method yields a very good 
approximation to the ideal of a Y-type neutral line on the source surface. 

In summary, our Implementation of the non-spherical source surface and least- 
squares method of imposing the coronal boundary condition there yield results that 
agree quite well with observation upon adjustment of only a few parameters. 
Further applications of the method to real and idealized solar conditions immediately 
suggest themselves and will be treated in future works of more specialized scope. 
We conclude, on the basis of results obtained so far, that the use of a norv 
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spherical source surface constitutes a significant procedural advancement in the 
modeling of coronal and heliospheric magnetic fields. 
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Appendix: Mathematical Detail 


We have deliberately kept the equations simple in the main text. For those 
readers who need or would prefer a mathematically explicit description of our 
procedures, we offer this Appendix. 

We assume that the scalar magnetic potential is given formally by equation (1) 
for Rq £ r ^ r^(9,^), but we recognize that the coefficients (g^, h^) and (g^, h^) 
are related as in equation (2) only for the special case of a spherical source 
surface. For the more general case of a non-spherical source surface we must 
differentiate a (cf. equation (6)) with respect to each of the g^i and each of the 
h^t in order to minimize a. We require specifically that 


da 


-ip' 


9fln 


f 3B dB 

®9n' ®9n' 


(Ala) 


and 

da f dB dB 

ah„, 8h„, ahn- 

o _ 1 

where dA=r |h 'r|~ d(cos9)d^ is the element of area on the source surface. 
Unless |n -rj is zero at some point on the source surface the area element dA is 
well defined. We intend to avoid such pathological cases, for both physical and 
mathematical reasons. We can make subsequent steps in the derivation of the 
(g^, h^) more compact by writing B=B + B , with B given by equation (7) and B 
given by 

N n ^ n 

B = -RqV 2 2 IrT] (gJ'cosm^ + h;f slnm^)p;f’(0), (A2) 

n=1m=0 ^ 


because we know that dB/dg™ =dB/dh™ =0. Thus, it follows from (A1) that 




dB _ dB 
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and 
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[BT^-(fiB)n.;^]dA 




&h'„, 


-X 


dB 


dB 


[(n • B) n • -r=fT - B • -Tjjjr] dA 
0h„"! ah„"«. 


(A3b) 


for each n' (=1, 2, 3 N) and m* (sQ, 1,2,.. ., n). We observe that 

n' 


9B ( f V m* 


(A4a) 


and 

®B / M"' m* 

^ = P™ («)»•"«.>]. (A4b) 

/.e., that both are independent of the expansion coefficients (g " , h^t ). Thus, it 

”2 ^ 

foiiows that (A3) constitutes a set of N +2N simultaneous linear algebraic 
equations for the (g™, h^). A more explicit representation of (A3a) is given in 

1 

equation (6). Equation (6) thus formally represents N (N + 3) separate linear 

equations for the N^+2N unknown parameters (a^> h^)- The other 'jN(N+1) 

required equations are obtained by replacing cos m'^ with sin m'^ (mVO) in 
equation (6) at each of the eight places where cos m'^ appears. 

The orthogonality properties of the harmonic functions reduce equation (6) and 
its complement to equation (2) for the case of a spherical source surface, /.e., for 
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tsRq and n^r. The more general and Interesting case o 1 a non-spherical source 
surface finds r=rg(5,^) and n^r In (6). The orthogonality properties of the 
harmonic functions are no longer useful, because the Integrand now contains vastly 
more complicated functions of 6 and ^ as well. However, we can express (6) and 
its complement In matrix form by ordering the coefficients (g^, h^) into an 

— p — 

(N + 2 N)-component vector C3 as follows: 

(C.J, Cg, C3, . . . , 

= (g^, g^J , h^ , ggt O2* ® 2 * ^ 2 ’ ^ 2 ’ ® 3 * • • • • (A6) 

Thus, priority In the sequence goes first to the smaller values of n over the larger, 
then to g over h, and finally to the smaller values of m over the larger. Similarly, 
the right-hand side of equation (6) and its complement, on which n" and m" are 
dummy variables only, can be ordered into an (N -t- 2 N)-component vector Kp by 
giving priority first to the smaller values of n' over the larger, then to equation (6) 
over its complement containing sin m'^ Instead of cos m'^, and finally to the smaller 
values of m* over the larger. With this ordering convention (or any other that is 
internally consistent) we can reduce equation (6) and its complement to the form 

N^+ 2 N 

s=1 

in which the elements of the A matrix are the integrals that appear as coefficients 
of g^ and h^ on the left-hand side of (6) and its complement. The matrix 
elements Ap^ must (of course) be arrayed in accordance with the convention 
chosen for ordering the vector components Kp and C^. By solving (A6) we obtain 
the set of expansion coefficients (g^, h^) in (1) that will minimize the mean-square 
tangential component of B over the source surface. This completes our basic 
objective. 
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It Is Important to ascertain whether the coefficients (g^, h^) thus determined 
are truly accurate for the case considered and reasonably stable against minor 
changes in N, in the (g", h^), or in the parameter a in equation (8), for example. 
Thus, it Is encouraging to know that our numerical code yields results in agreement 
with equation (2) for a=0, in agreement with expectations dictated by symmetry 
(see Section 3) in the case of a photospheric dipole for general a, and in 
agreement (see Section 4) with the numerical results of Schulz et al. (1978) for 
the special case of a photospheric dipole with a=1. When we obtain the 
coefficients (g^, h^) from real magnetic data rather than from an idealized model 
of the photospheric B field, then (of course) there is no standard set of 
coefficients (g^, h^) against which to compare those calculated by minimizing <r 
over a non-spherical source surface. Indeed, standard sets of coefficients 

(g^, h") are available only for the simplest of idealized models (Schulz et a/., 
1 978). However the stability of coefficients (g^> h^) determined numerically by 
minimizing ff, either in a realistic case or In an idealization more complex than the 
photospheric dipole case, can be Judged by an inspection of the A matrix that 
appears In (A6) and whose elements are the Integrals that appear as coefficients 
of the unknown quantities (g", h^) on the left-hand side of equation (6) and its 
complement. For the (g^, h^) determined by solving these linear equations to be 
stable against small changes in r^(9,^) or N, the matrix (A) of the coefficients must 
be well-conditioned, i.e„ free of eigenvalues that either vanish or are too small 
because of a failure (or near^failure) of the rows of A to be linearly independent of 
each other. Our experience suggests that we have little to fear In this regard. 
While the A matrix is strictly diagonal only in the limit of a spherical source surface 
(a=0), it Is nearly diagonal In the other cases that we have examined. In other 
words, the elements Ap^ In (A6) become progressively and substantially smaller in 
magnitude as one steps away from the main diagonal (p=s). This property is 
generally sufficient to assure that a matrix is well-conditioned. Moreover, we have 
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encountered no obvious difficulties in seeking solutions of (A6). The solutions for 
the (g™> h™) seem to vary continuously with a, for example (see Figure 3), and 
weakly with N. Thus, we find overall that our method is well-behaved, and we have 
full confidence in it. 

Having thus obtained the optimal coefficients (g^, h^) for use in (1), we can 
proceed to plot field lines inside the source surface, /.e., for r ^ This 

IS easy enough to do, starting either from a footpoint at r=Ro or from a point on 
the source surface, because the trajectory of a field line is defined locally by the 
direction of B. Having obtained the coefficients (g^, h^) and (g 
calculate the components (B^, B^, B^) of B=-VV from (3). As usual the field line 
must satisfy the equations 



where B=|B| and ds is the element of arc length along the trajectory. We use a 
Runge-Kutta method lo be certain that the computed trajectory truly tracks the 
field line (despite the finite step size). Outside the source surface we plot field 
lines with the aid of (13), thus imposing a boundary condition that B is exactly 
normal to the source surface. In other words, we do got “continue" the direction 
of the numerically computed magnetic field across the source surface from inside to 
outside, because there might be a significant error locally in doing so (despite the 
global minimization of <r). Thus, the footpoint of the interplanetary segment of a 
field line coincides with the point at which the coronal segment of the same field 
line touches the source surface from the inside. We accept the small discontinuity 
in h X B that is inevitable in such a construction when the source surface is non- 
spherical, but we would assume the continuity of h * B across the source surface in 
order to calculate the magnitude of |B| in interplanetary space (Schulz et a/., 
1978). The actual calculation of |B| along a field line outside the source surface 
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involves numerical application of differential geometry and Is not yet completed. 
Thus, we have been content In the present work to obtain the direction of B from 
the fact that 


B 


(* X + y y + ± £) 


(A8) 


and to trace the field lines themselves by means of equation (13). Field lines 
would be rectilinear outside the source surface except for the garden-hose effect 
caused by the rotation of the sun (Parker, 1968). indeed, at heliocentric 
distances rj^lORQ the field lines described by equation (13) deviate negligibly 
from rectilinear projections of the outward normals n from the various points on the 
source surface, and we have used tMs fact to simplify the computations underlying 
Figures 3, 7, and 8. However, the spiral configuration of B is very important at 
interplanetary distances, as Is well known. Thus, ft has been absolutely essential 
to use equation (13) in its full generality for the purpose of locating the 
Interplanetary current sheets In Figure 1 2 (especially so In Figure 1 2b). 


52 



References 


Adams, J. and Pneuman, G.W.: 1976, Solar Phys., 46, 185. 

Altschuler, M.D. and Newkirk, G., Jr.: 1969, Solar Phys., 9, 131. 

Altschuler, M.D., Levine, R.H., Stix, M. and Harvey, J.W.: 1976, Solar Phys., 61, 346. 

Bohlin, J.D.: 1977, in J.B. Zirker (ed.), Coronal Holes and High Speed Solar Wind 
Streams, Colo. Assoc. Univ. Press, Boulder, pp. 27-69. 

Burlaga, L.F., Behannon, K.H., Hansen, S.F. and Pneuman, G.W.: 1978, J. Geophys. 
Res., 83, 41 77. 

Howard, R.: 1977, Ann. Rev. Astron. Astrophys., 15, 153. 

Hundhausen, A.: 1 970, Rev. Geophys. Space Phys., 8, 729. 

Jackson, B.V. and Levine, R.H.: 1 980, Solar Phys., in press. 

King, J.H. (ed.): 1975, Interplanetary Magnetic Field Data, World Data Center A for 
Solar-Terrestrial Physics Report UAG-46, NOAA, Boulder, Colo. 

Levine. R.H.: 1977, Astrophys. J., 218, 291. 

Levine, R.H.: 1978, J. Geophys. Res., 83, 4193. 

Levine, R.H. and Altschuler, M.D.: 1974, Solar Phys., 36, 346. 

Levine, R.H., Altschuler, M.D. and Harvey, J.W.: 1977, J. Geophys. Res., 82, 1061. 

Munro, R. and Jackson, B.V.: 1977, Astrophys. J., 213, 874. 

Nerney, S.F. and Suess, S.T.: 1975, Astrophys. J., 196, 837. 

Parker, E.N.: 1958, Astrophys. J., 128, 664. 

Pneuman, G.W. and Kopp, R.A.: 1971, Solar Phys., 18, 258. 

Poletto, G., Vaiana, G.S., Zombeck, M.V., Krieger, A.S., and Timothy, A.F.: 1975, Solar 
Phys., 44, 83. 

Riesebieter, W. and Neubauer, F.M.: 1979, Solar Phys., 63, 127. 

Schatten, K.H.; 1968, Nature, 220, 1211. 

Schatten, K.H., Wilcox, J.M. and Ness, N.F.: 1 969, Solar Phys., 9, 442. 


53 



Schulz, M.: 1073, Astrophys. Space Sci., 24, 371. 


Schulz, M., Frazier, E.N. and Boucher, D.J., Jr.: 1978, Solar Rhys., 60, 83. 

Smerd, S.F. and Dulk, G.A.: 1971, in R. Howard (ed.), IMJ Symposium No. 43, Solar 
Magnetic Fields, Reidel, Dordrecht-HoHand, pp. 618-641. 

Smith, S. and Schatten, K.H.: 1 970, Nature, 226, 1 1 30. 

Sturrock, P.A. and Smith, S.M.: 1 988, Sdar Rhys., 6, 87. 

Suess, S.T., Richter, A.K., WInge, C.R. and Nemey, S.F.: 1077, Astrophys. J., 217, 
296. 

Svalgaard, L, Duvall, T.L., Jr., and Scherrer, P.H.: 1 978, Solar Rhys., 56, 226 . 
Svestka, Z., Solodyna, C.V., Howard, R. and Levine, R.H.: 1977, Solar Rhys., 65, 359. 
Uchida, Y., Altschuler, M.D. and Newkirk, G.A., Jr.: 1 973, Solar Rhys., 28, 496. 


54 



End of Document 



